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ABSTRACT 

Compact object clusters are likely to exist in the centre of some galaxies because of mass seg- 
regation. The high densities and velocities reached in them deserves a better understanding. 
The formation of binaries and their subsequent merging by gravitational radiation emission 
is important to the evolution of such clusters. We address the evolution of such a system in 
a relativistic regime. The recurrent mergers at high velocities create an object with a mass 
much larger than the average. For this aim we modified the direct NBODY6++ code to include 
post-Newtonian effects to the force during two-body encounters. We adjusted the equations 
of motion to include for the first time the effects of both periastron shift and energy loss by 
emission of gravitational waves and so to study the eventual decay and merger of radiating 
binaries. The method employed allows us to give here an accurate post-Newtonian descrip- 
tion of the formation of a run-away compact object by successive mergers with surrounding 
particles, as well as the distribution of characteristic eccentricities in the events. This study 
should be envisaged as a first step towards a detailed, accurate study of possible gravitational 
waves sources thanks to the combination of the direct Nbody numerical tool with the imple- 
mentation of post-Newtonian terms on it. 
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1 INTRODUCTION 

It is nowadays well established that most, if not all, galaxies har- 
bour a supermassive black hole (SMBH) in their centre with a mass 
of some 10 6 ~ 9 M(^ (see e.g . the recent review b ylFerrar ese & Forcj 
l2005t iFerrarese et aljEool iKormendy & Gebhardtl|200lh . There 
are also signs for masses of 10 6 Mm iGreene^^Hq|200a) . In the 
case of our Galaxy this is even i mperative; an SMBH with a 
mass of about ~ 3 - 4 x W 6 M p, <Eckart et all2002tlGh"ez et alj 
|2000L 120031 ISchodel et al]|2002h must be ensconced in its cen- 
tre. If one extends the correlation between the SMBH mass and 
the stellar velocity dispersion of the bulge of the host galaxy (the 
Mm — a correlation) observed for galactic nuclei iGebhardt et alj 
1200(1 IFerrarese & MerritfeOOd) to smaller systems, like globular 
clusters, one should expect intermediate mass black holes (IMBH) 
with masses of between 10 3 — 1O 4 M0 to be lurking in the cen- 
tres of such stellar clusters. There are observations of M15 in the 
Milky Way or Gl in M31 iGerssen et alj|2002l : IGebhardt et all 
2002; van der Marel et alj|2003l) which are compatible with this 
possibility, but N— body models of these clusters have been made 
which do not require the presence of an IMBH iBaumgardt et alj 
12003d) . 

The densities observed in the central region of galaxies, where 
these very massive objects are located, are very high and may even 
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exceed the core density of globular clusters by a factor hundred 
(about 10 7 - 10 8 M o pc~ 3 for the Galactic Centre, for instance) 
and thus make them very special laboratories for stellar dynamics. 

On the other hand, it is not strictly excluded that the central 
mass concentrations are not massive black holes (MBHs). Mass 
segregation creates a flow of compact objects like ne utron stars 
or stel lar black holes to the central parts of the cluster <Leell987t 
iMiralda-Escude & Goul ol200d) and may constitute there a cluster. 
This could mimic the effect of the MBH, and thus give an alterna- 
tive explanation of the pro perties of clusters tha t have gone core- 
collapse like M15 and Gl iGebhardt et all2002HBaumgardt et alj 
l2003alrilvan der Marel et al.l2003l) . On the other hand, MB Hs are 
favoured in the cas e of galaxies, in particular the Milky Way iMaozl 
ll998tlMillej2005h . 

For the case of a globular cluster it has been studied that 
stellar black holes are probably ejected from the system. Stellar 
black holes should form three body binaries and kick each other 
out of the cluster jPhirmex& Si^urdssorj LL99lt |Kulkarru_et_alJ 
1993 ; ISigurdsson & Hernauistll993UPortegies Zwart & McMillan! 
2000). Nonetheless, if the velocity dispersion is high enough, then 
binaries will not be created due to three body encounters, as in the 
classical case considered before, but to gravitational waves emis- 
sion during two-body encounters. A simple way to understand this 
is that the components of a binary merge before a third particle or 
a second binary comes in sufficiently close to interact with them so 
as to eject the binary or one of its. Thus, ejections cannot happen 
in such a scenario. As a matter of fact, for velocity dispersions of 
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> 300 km/s the merging time in clusters with two mass compo- 
nents is already shorter than the required time between interactions 
before a th ird particle or a second binary can bring about an ejection 
fcxell995h . 

Relativistic stellar dynamics is of paramount importance for 
the study of a number of subjects. For instance if we want to have 
a better understanding of what the constraints on alternatives to su- 
permassive black holes are; in order to canvass the possibility of 
ruling out stellar clusters, one must do detailed analysis of the dy- 
namics of rela tivistic cluste rs and determine in particular the core 
collapse time jMillej2005l) . Also, when we want to more compe- 
tently dive into the formation of MBHs, learn how the dynamics 
around them is, for instance to estimate captures of compact ob- 
jects on a central SMBH via extreme mass ratio inspiraling, or pe- 
ruse a system of many supermassive black holes etc the inclusion 
of relativistic effects is constitutive. 

Our current work includes the study of stars on relativistic or- 
bits around a SMBH, so as to be able to estimate captures of com- 
pact objects on a central SMBH via extreme mass ratio inspiraling 
and binary evolution of two SMBHs. 

Efforts to understand the dynamical evolution of a stellar clus- 
ter in w hich relativ i stic effects may be important have been already 
done b v lLe3h987t).lQuinlan & Shanir3<1989h . lOu'inlan & Shaplrol 
1 19901) and lLeeH 19931) . In his work. lLeeU 19931) (MHL93 from now 



onwards) addressed the problem of the dynamical evolution of a 
cluster composed of compact objects by, with some approxima- 
tions, adding an estimate of the gravitational wave emission term 
correction to NBODY5 (see section 3). Nevertheless, he neglected 
the 1P M and 2PM t erms and made use of the formalism intro- 
duced bv lPetersliT964l) . possibly because of computational bourns. 
The computation of the PM corrections is CPU-consuming, for we 
have to compute both, the accelerations and their time-derivatives 
(see next section). Also, NBODY5 is not suitable for supercomput- 
ers or special purpose GRA PE hardware; here either NBODY6+ + 
or NBODY4 have to be used lSpurzemll999l : IS. J. Aarsethll99^) . 

In this work we describe a new tool that allows us to address 
this problem in a much more rigurous way than done in the existing 
literature, including deviations from th e Newtonian fo rmalism of 
the standard direct NBODY6+ + code |S purzem|| l999h . based on 
Aarseth' s direct NBODY codes IS. J. Aarsethll999h . We modified it 
in order to allow for post-Newtonian (PM) effects, implementing 
in it the 1PM, 2PM and 2.5PM corrections without any further 
approximation th an those indw elling to the calculation of the VM 
terms themselves JSoffelll989h . 

In Section 2 we give a brief description of the method and 
of the implementation of the PM terms into a standard NBODY 
code. An analysis of the formation and evolution of a particle that 
gains more and more mass from successive mergers in the system 
(the "runaway particle") is made in Section 3 and, to conclude, in 
Section 4 we make a summary and discussion of the main results 
obtained. 



2 METHOD: DIRECT SUMMATION NB ODY WITH 
POST-NEWTONIAN CORRECTIONS 

The version of direct summation NBODY method we employed 
for the calculations, NBODY6++, includes the KS regularisa- 
tion. This means that when two particles are tightly bound to 
each other or the separation among them becomes smaller dur- 
ing a hyperbolic encounter, the couple becomes a candidate for 
a in order to avoid problematical small individual time steps 



iKustaanheimo and Stiefel 119651) . We modified this scheme to al- 
low for relativistic corrections to the Newtonian forces by expand- 
ing t he acceleration in a series of powers of 1/c in the following 
way iDamour & Dereullell98il ; ISoffelll989l) : 
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where a is the acceleration of particle 1, a — —Gm2n/r 2 is 
the Newtonian acceleration, G is the gravitation constant, rti\ and 
m 2 are the masses of the two particles, r is the distance of the 
particles, n is the unit vector pointing from particle 2 to particle 
1, and the 1PM, 2PM and 2.5PM are post-Newtonian correc- 
tions to the Newtonian acceleration, responsible for the pericen- 
ter shift (1PM, 2PM) and the quadrupole gravitational radiation 
(2.5PM), correspondingly, as shown in Eq. Q. The expressions 
for the accelerations are: 
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In the last expressions and v 2 are the velocities of the parti- 
cles. For simplification, we have denoted the vector product of two 
vectors, x_ x and x 2 , like X\X2- The basis of direct NBODY4 and 
NBODY6++ codes relies on an improved Herm it integrator scheme 
iMakino & AarsetlJI 19921 IS. J. Aarsethll 19991) for which we need 
not only the accelerations but also their time derivative. These 
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derivatives are not included in these pages for succinctness. We in- 
tegrated our correcting terms into the KS regularisation scheme as 
perturbations, similarly to what is done to account for passing stars 
influencing a KS pair. Note that formally the perturbation force in 
the KS for malism does no t need to be small compared to the two- 
body force <Mikkolall997t) . If the internal KS time step is properly 
adjusted, the method will work even for relativistic terms becoming 
comparable to the Newtonian force component. 



3 DYNAMICAL EVOLUTION OF A CLUSTER OF 
COMPACT OBJECTS 

3.1 The initial system and units 

The units used in our models correspond to the so-called TV-body 
unit system, in which G = 1, the total initial ma ss of the stel- 
lar cluster is 1 and its in itial total energy is -1/2 <Henonlll97lt 
iHeggie & Mathieull98d) . The sy stem was ch osen to be initially to 
be indentical to that employed bv lLee|jl993l) :i.e. a spherical clus- 
ter with a number of compact stars A/* = 10 3 of identical mass 
m. These were distributed in an isotropic Plummer sphere, which 
means that the phase-space distribution function is proportional to 
\E\ 7 / 2 , where E is the energy per unit mass of one star. The den- 



sity profile is thus p(r) = po(l + (r/Rp\) 



5/2 



where Rpi is 



the Plummer scaling length. For such a model the TV-body length 
unit is U\ = 16/(37r) 7?pi. 

In the situations considered here, the evolution of the cluster 
is driven by 2-body relaxation. A natural time sca le is the (initia l) 
half-mass relaxation time. We use the definition of lSpitzeJll987l) . 
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For instance, for a Plummer model, the half-mass radius is R1/2 = 
0.769 Wi = 1.305 Rpi. M c i is the total stellar mass and In A = 
In (yN) is the Coulomb logarithm. 

For the situation considered in this work, the square ratio of 
the central velocity dispersion <r cen to the speed of light c, 
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is big enough, so that we can expect that relativistic effects play 
a noticeable role in the evolution of the system. For this aim, we 
chose CT con to be ~ 4300 km/s. G is the gravitational constant, 
7l cl is the radius of the cluster and Tlf chv , = 2GM cl /c 2 is the 
Schwarzschild radius of the cluster. 

In our calculations the VAf terms are acting all the time dur- 
ing the calculations but obviously become important only when ve- 
locities are high. Our criterion for particles to merge is that they 
reach their common Schwarzschild radii 7?.g c hw; i.e. the sum of 
their Schwarzschild radii. This is of course approximative because 
the VM treatment breaks down when particles are that close (and 
v ~ c), but this should not matter, for the merging phase is much 
faster than any stellar dynamical time. The gravitational recoil, the 
expected lose linear momentum in asymmetric systems in which 
the merger remnant receives a kick from the gravitational waves 
emission obviously does not show up in our models, because we 
truncate the series at 0(cT b ) and it is only to be treated as an effect 
of higher-order terms. 



3.2 Formation of a run-away body 

Even though we started with a single-mass stellar system, the 
masses of some objects in the cluster increased by relativistic merg- 
ers. In Fig. Q we survey the time evolution of the mass increase. 
We find a number of mergers that lead to the variation of the initial 
single mass situation. The particle masses increase after the rel- 
ativistic merging events, since we are assuming that the particles 
merge perfectly when they reach the distance of their 7£schw (see 
above). We find the formation of a runaway particle that reaches 
almost six percent of the initial total mass by the end of the simu- 
lation, see Fig. 0. We denoted the mass of runaway body by red 
crosses and the mass of other mergers by blue crosses. 

One can observe that the runaway body dominates the system 
after its fast growing phase around 300 time units, which is ap- 
proximately the moment at which the core collapse of the system 
happens, as we can see in Fig. Only some merger events which 
are independent from the runaway body can occur after this phase. 
This fast growing phase occurs at the core collapse of the system 
iMevlan & Heggiell997l) . In Fig. 13 we follow the evolution of the 
so-called Lagrangian radii of the system, spheres containing 1%, 
5%, 10%, 20%, 30%, 50%, 70%, 90% and 100% of the total mass 
of the cluster; the centre of the cluster is defined to be the centre 
of the mass density. Since the runaway particle is included, and in 
the end its mass reaches 5% of the total initial mass of the cluster, 
the curves corresponding to 1% and 5% roughly correspond to its 
evolution. We observe that the runaway stops the core collapse and 
allows for a expansion. 

The process of mergers translates directely into a production 
of energy in the central regions of the cluster. The centre adapts 
to supply the cluster with the same amount of energy that it can 
obtain via relaxation, and this amount is determined by the large- 
scale structure. 

According to, for instance, the table given in lFreitag & Benz| 
1200 ll) . the standard value for the core collapse time is of roughly 
~ 15 — 20 times the half-mass relaxation time T r h- We find nev- 
ertheless that the core collapse time is t c c ~ UTa, with a valu e 
of 7 = 0.11 in the Coulomb logarithm iGiersz & Heggiell 19941) . 
which clearly suggests that the VM terms accelerate the collapse. 
This can be seen more clearly in Fig. J2j> which corresponds to the 
same simulation but without making use of relativistic corrections. 
There we can see that t cc ~ 380 ~ 14T r h. 

In Fig. {3j we show the evolution of the the runway particle 
mass normal ised to the mass contained in the core of the cluster, 
defined as in lCasertano. S. & Hut. Pl ll985t). The mass of the run- 
way particle can grow only up to the core mass. The core mass 
continuously decreases as the core collapse proceeds. We see this 
in the figure, where the runaway particle grows and saturates to the 
core mass after ~ 1200 time units. 

The evolution and formation of the runaway particle mass is 
not as fast as it was in MHL93, as we can see in his Figure 5. For 
our simulation the sudden jump in the growth of the mass comes 
in slightly later and is smoother, reaching final values for the run- 
away particle mass of about three times smaller than in MHL93. 
The reasons for the differences are to be attributed to the following: 
MHL93 calculated the influence of the 2.5V Af term on the orbits in 
an unperturbe d pair and ma de them merge after a decay timescale, 
following the Peters ( 1964) formalism. This requires the assump- 
tion that particles move along their orbits on an ellipsis, only valid 
when they are very far from the relativistic regime. On the other 
hand, we implemented the 2.5V M term in the code itself, so that 
the relativistic corrections are a natural feature whose influence on 
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Figure 1. Time evolution of merging masses. The formation of the runaway 
particle is about the time of the cluster core collapse. Fore more details see 
text 
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Figure 2. Evolution of the Lagrangian radii corresponding from the bottom 
to the top to 1%, 5%, 10%, 20%, 30%, 50%, 70%, 90% and 100% of the 
total mass 



the evolution of the system comes in when the velocities of the 
stars become high enough. The influence of the IV N and 2VJ\f 
terms, corresponds to the conservative phase evolution of the orbit 
and cannot be relevant because they do not change its energy and 
angular momentum. 
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Figure 3. Evolution of the runaway particle mass in units of the core mass 
(at the same time) 



1 VN, 2 VN terms modify the extrinsic features of the orbits (e.g. 
the orientation) but do not affect their intrinsic parameters (like fre- 
quency), we therefore can expect their effect to be averaged out 
during the evolution of the system and not influence the mergers 
rates. One should thus attribute the differences to the approach he 
made, somehow inadequate for the velocity regime considered. 

This study should be envisaged as successful test test of the 
code, which showed to be robust. This tool can be applied to other 
astrophysical scenarios which require a post-Newtonian treatment. 
This includes on-going work, as e.g. captures of compact objects by 
a supermassive black hole in a galactic centre, also known as ex- 
treme mass ratio inspirals (EMRIs). One of the fundamental aims 
is to rigurously explore the parameter space, so that we can provide 
the LISA data analysis community with realistic estimates of, for 
instance, the eccentricity, mass ratio etc at the beginning of the fi- 
nal merger, when the smaller compact object enters the LISA band. 
An assumption for the initial parameter space is necessary in or- 
der to develop waveform "banks" for this kind of events. One must 
note here that the inclusion of the 1 VN and 2 VN terms is very 
relevant, for ressonant relaxation (or Kozai) effects, which may in- 
crease the rate of inspiral significantly, may be strongly affected by 
by relativist ic precession and thus have an impact on the number 
of captures iHopman & Alexanderll200rj|Kozailll962l) . The inclu- 
sion of higher-order VN terms is also part of current study and will 
also shed light on other aspects of this subject (spin-spin coupling, 
spin-orbit interaction and radiation recoil). 



4 CONCLUSIONS 

In this work we have presented a study of the formation and evo- 
lution of a runaway particle in a dense cluster of compact ob- 
jects -which initially had the same mass- as a result of relativistic 
mergers. We employed a modified version of the direct summation 
NBODY6 code in which we have implemented the 1VN, 2VN, 
and 2.5 VN terms to take into account post-Newtonian corrections 
to the standard NBODY Newtonian acceleration. 

The runaway particle reaches in the end of our simulations 
~ 6% of the initial total stellar mass of the cluster. We have also 
compared our work to a previous result ba sed on a more approx- 
imative scheme, the approach described in lPeters) J 19641) and we 
have found out that the net result is that the growth of the run- 
away particle in the study of MHL93 is ~ 3 times larger. Since the 
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